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SECTION  I 
INTRODUCTION 

The  calculation  of  nuclear  fallout  doses  and  resulting  casualties  has  been 
practiced  by  many  people  for  many  years  using  many  different  algorithms.  One 
of  the  difficult  areas  in  fallout  assessments  is  the  specification  of  wind  fields 
required  to  calculate  the  paths  of  fallout  particles.  Many  studies,  particularly 
those  with  large  numbers  of  weapons  and  targets,  use  mean  winds  or  some  "effective" 
wind.  This  is  to  allow  calculation  of  fallout  effects  of  large  scenarios 
within  a reasonable  amount  of  time.  Other  studies,  with  much  more  modest 
scenarios,  have  used  tabulated  winds  to  determine  fallout  effects.  Large  data 
bases  have  been  built  containing  measured  wind  fields  over  the  Continental 
United  States.  The  present  effort  presents  a new  methodology  to  arrive  at  wind 
fields  for  modest  scenarios  and  takes  a brief  look  at  the  adequacy  of  mean  wind 
methods  as  compared  to  using  more  realistic  wind  fields. 
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SECTION  II 

THE  GENERATION  OF  THREE-DIMENSIONAL  CORRELATED  WINDS 


Let  the  wind  be  represented  at  a point  by  a complex  number 


w(x,y,z)  = wxlx,y,z)  + i wy(x,y,z) 


(1) 


where 


x = west  to  east  coordinate 
y = south  to  north  coordinate 
z = vertical  coordinate 


w = west  to  east  wind 
x 

wy  = south  to  north  wind 


It  is  assumed  that  the  mean  wind  components  are  not  contained  in  equation  (1), 


hence,  w and  w are  zero  mean  random  variables.  After  w has  been  generated, 
x y 

the  mean  components  can  be  added  in  to  get  total  winds.  It  is  furthermore 


assumed  that  w and  w gaussian  distributed  random  variables  and  the  second 
X V 


order  statistical  processes  with  respect  to  w and  w are  stationary  and  homo- 

x y 


geneous.  With  these  assumptions,  the  statistics  of  the  winds  are  completely 
described  by  the  complex  correlation  function 


G(p,n,e)  * w*(x,y,z)  w(x  + p,  y + n,  z + r) 
and  the  asymmetric  correlation  function 


(2) 


F(p,n,e)  3 w( x,y ,z)  w(x  + p,  y + n,  z + e)  (3) 

Let  us  represent  a wind  pattern  as  a three-dimensional  fourier  integral. 


where 


w(x,y'2)  * itf  { dk*  c(k*vt*)e  1<v  * v * k‘2>  <4) 

c(kx.ky,kz)  dx  j dv  J dz  W^x,y,z^e  "^kxX  + kyy  + kzz> 


(5) 


! . J r 
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It  can  be  shown  that 

c*(kx,ky,kz)  c(kx,ky,kzJ  = 

J'dp  J“dr\  J'de  G(p,n,e)  e ^kx°  kyn 

c(kx,ky,kz)  c(-kx,-ky,-kz)"  = 

J'dp  Jdr)  J'de  F(p,n,e)  e ^kxp  + kyn 

and  that  all  other  combinations  of  k ,k  , and  k are  zero.  Equations  (6)  and 

x y l 

(7)  show  that  fourier  coefficients  correlate  only  in  pairs  mirrored  through  the 
origin  of  fourier  space.  It  can  also  be  shown  that  if  w(x,y,z)  is  gaussian 
distributed,  then  so  are  the  fourier  coefficients.  Finally,  it  follows  that  if 
two  fourier  coefficients  are  gaussian  distributed  and  uncorrelated,  then  they 
are  also  statistically  independent.  By  going  to  fourier  space,  the  sampling 
problem  has  been  reduced  to  sampling  pairs  of  statistically  independent  complex 
random  numbers. 

In  practice,  the  fourier  integral  expressions  are  approximated  by  discrete 
transforms  evaluated  using  fast  fourier  transform  algorithms.  Equation  (4)  is 

now  written  as  / , 

ijil*  + hi.  + illA 

w(x,y,z)  = a-.-  f .■■■  EE  E cikJ  e \x  LV  Lz  / (8) 

01 xyLz  j k i Jia 

where  -Lx  < x < Lx,  -L  < y < Ly,  and  -Lz  < z < Lz>  Equations  (6)  and  (7)  are 
now 


+ kze)  (6) 

+ kze)  (7) 
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The  only  basic  numerical  restrictions  on  using  equations  (8)-(10)  is  that  the 
spatial  resolution  adequately  resolve  G(p,n,e)  and  F(p,n,e)  and  that  the  afore- 
mentioned functions  become  small  at  the  spatial  mesh  boundaries. 

Let  us  now  look  at  a specific  pair  of  coefficients,  and 

Let 


cjki ■*."  bl 

(11) 

C-j-h-t  * a2  * ' b2 

(12) 

C*jkn  Cjkn  = K 

(13) 

CV  ■ L * 1 H 

(14) 

From  equations  (11)-(14)  with  j,  k and  s.  not  equal  to  zero 

ai2  = b!2  = a22  = b22  = K/2  (15) 


a1b1  - a2b2  = 0 


(16) 


aia2  = L/2 


ajb2  = a2bj  * M/2 

The  joint  probability  distribution  function  of  , a2 , b^,  and  b2  is 

r 


(17) 

(18) 


P(ai , bi , a2,  b2 ) = C exp  / 


“ 

1 

- 

-K  0 -L  -M 

*1 

1*1  bi  a2  b2| 

0 K -M  L 

bl 

-L  -M  K 0 

a2 

M L 0 K 

b2 

(K2  - L2  - M2) 


j 


v. 
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If  a new  set  of  variables  is  defined  such  that  if  yfi2  + M2  = 

-M/V 


and 


a. 

l 

b, 

l 

a. 

= 

2 

b 

2 

1 0 

0 1 

L/V  M/v^ 

M/\/  -L/V 


■L/V” 

-M/V 

1 

0 


L/V 

0 

1 


b2 


1 

'2(ap  + bp)  2(a’2  + bp)' 

+ 

‘ 2 

. K + V K - V . 

(20) 


then  after  considerable  matrix  manipulation  it  can  be  shown  that  the  joint 
probability  distribution  of  the  new  variable  set  is 


P'  (a- , b;,  a',  bp  = C'  exp 


(21) 


The  new  variable  set  can  be  easily  sampled  from  a gaussian  distribution  because, 
as  shown  in  equation  (21),  they  are  independent.  After  the  new  variable  set  has 
been  sampled,  the  original  set  can  be  retrieved  using  equation  (2)  and  Cjkjl  and 
can  be  calculated  from  equations  (11)  and  (12).  For  j = k = l = 0,  the 
sampling  procedure  is  slightly  different  because  aT  = a2  and  b:  = b2 . As  before, 
a new  set  of  independent  coefficients  can  be  found  that  can  be  sampled  indepen- 
dently, Let 


a,  = 


(K  + L)/2]X 


/ 2 


(a{  + b{) 


bi  [(K  - D/z]1/2  (a{  + b{) 

Then  the  joint  probability  distribution  for  aj  and  b{  is 


P(a{ , b{ ) = C exp 


J 


2\ K2  - L2  - M(K2  - L2) 
K2  - L2  - M2 


l/2) 


K2  - L2  + M(K2  - L2) 
K2  - L2  - M2 


1/2 


(22) 

(23) 


(24) 


Sample  values  for  aj  and  b1  can  be  found  by  sampling  aj  and  b(  according  to 
equation  (24)  and  then  evaluating  equations  (22)  and  (23). 

The  procedure  for  generating  a sample  zero  mean  wind  pattern  can  be 
summarized  as  follows.  The  Cjki  for  all  j,  k,  and  l are  sampled  using  the 
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procedures  described  above.  The  wind  field  is  then  found  by  evaluating  equation 
(8).  For  the  total  wind  fiel^,  the  mean  winds  can  be  added  to  w(x,y,z). 

It  was  earlier  assumed  that  the  second  order  statistics  of  the  wind  field 
were  homogeneous.  This  restriction  can  be  easily  relieved,  particularly  with 
respect  to  the  standard  deviation  of  the  wind  components.  In  practice  the  zero 
mean  wind  field  is  generated  assuming  unit  standard  deviation  for  both  wind 
components.  The  final  winds  are  calculated  by  multiplying  the  randomly  generated 
wind  components  by  their  local  standard  deviations.  Other  inhomogeneties  can  be 
approximated  in  a similar  manner  if  required. 
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SECTION  III 

THREE-DIMENSIONAL  WIND  CORRELATION  FUNCTIONS  FOR  GEOSTROPHIC  WINDS 


The  wind  modeling  used  in  this  report  requires  the  wind  correlation  functions 
Buell  (refs.  1 and  2)  has  developed  several  candidate  horizontal  wind  correlation 
functions  which  do  a very  reasonable  job  in  matching  available  measurements. 

There  are  two  functions  required  to  specify  the  wind  correlations.  The  first  is 
the  longitudinal  correlation  function  which  describes  the  wind  correlation  at 
two  points  where  the  components  sampled  are  parallel  to  the  vector  connecting 
the  two  points.  The  second  function  is  the  transverse  correlation  function  where 
the  sampled  components  are  perpendicular  to  the  connecting  vector.  The  functions 
used,  assuming  unit  standard  deviation  are 

ru,(r)  =e'a2r2/2 

(25) 

ru,(r)  = O ' a2r2)e-a2r2/2 

where 

r = (p2  + n2 ) !/2 

(26) 

a = correlation  length 

The  correlation  length,  a,  is  taken  as  1.2  x 10*3  km"1. 

Equations  (1)  and  (2)  are  easily  converted  to  Cartesian  coordinates 
F(p,n)  and  G(p,n)  functions 

and  then  to 

R(p»n)  = e a r ^ [2  - a2  (p2  + n2)] 

(27) 

S(p,n)  = e x r [at2(n2  - p2)  + 2i  a2p.J 

(28) 

For  the  vertical  correlation,  an  exponential  correlation  function  is  used. 

T(e)  = e_Bc  (29) 

This  form  is  supported  by  Kochanski  (ref.  3)  who  found  that  B varied  from  0.08  to 
0.14  km"1  during  the  Fall,  Winter  and  Spring  over  the  Continental  United  States. 
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During  the  sunnier  B varied  from  0.11  to  0.15  km*1.  For  present  purposes,  B was 
taken  to  be  0.10  km*1 . 

The  total  correlation  functions  are 


G(p,n>e)  = R(p,n)  T (e) 


(30) 


F(p  ,n  ,e)  = S(p ,n)  T (e)  (31 ) 

The  discussion  on  generating  three-dimensional  wind  patterns  showed  that  the 
fourier  transforms  of  equations  (30)  and  (31)  are  the  quantities  that  are  used 
directly.  The  transforms  are 


/ dp  J dn  J 


dp  G( 


P >n  »e 


)e 


-i(kxp+kyn+kze) 


-fk  2+k v2)/2a2 
4-rrBe  v x y ' 

a4  .V  + k 2^ 

V z / 


(32) 


f ^P  y*dn  J 


de  F(p,n,e)e 


-1(Vky+kzc) 


4TrBe  v y 
a4  (B2  + kz2) 


kx2  ‘ 21  kx  ky) 
(33) 


The  first  and  second  order  statistics  for  the  winds  used  in  this  study  are 
from  the  Joint  Strategic  Target  Planning  Staff  Climatic  Wind  Data  Base. 

This  data  base  contains  the  mean  west-east  winds,  the  mean  south-north  winds, 
the  major  and  minor  axis  of  the  elliptical  standard  deviation,  and  the  rotation 
angle  of  the  ellipse  measured  counterclockwise  from  the  west-east  direction  for 
points  on  a latitude-longitude  array  and  for  altitudes  ranging  from  sea  level  to 
about  20  km.  The  standard  deviation  ellipses  over  the  U.S.  are  near  circular 
with  the  major  axis  averaging  1.2  times  the  minor  axis.  The  average  of  the 
rotation  angle  is  near  zero  and  the  standard  deviation  is  about  25°.  Because  of 
these  properties  of  the  standard  deviation  ellipses,  the  rotation  angle  is  ignored 
except  to  calculate  the  west-east  and  south-north  wind  standard  deviations  used 
to  normalize  the  random  wind  mesh. 
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SECTION  IV 

FALLOUT  TRANSPORT  MODEL 

There  are  four  basic  elements  to  any  transport  model.  They  are  the  calcula- 
tion of  the  amount  and  time  dependence  of  the  activity  created,  the  particle 
size  distribution,  the  activity  per  particle  of  a given  size,  and  the  fallout 
cloud  dynamics.  The  first  three  are  modeled  with  simple  models  that  are  physi- 
cally reasonable  but  are  not  represented  as  particularly  rigorous.  The  early 
time  cloud  dynamics  is  also  modeled  relatively  simple.  It  is  the  late  tine 
transport  where  the  natural  wind  field  dominates  that  is  emphasized.  Each  of 
the  basic  elements  are  discussed  in  turn. 

The  total  amount  of  activity  created  is  assumed  to  be  1500  roentgen  per  hour 
per  kT  per  square  mile  at  1 hour  after  burst.  If  the  burst  is  at  the  ground 
level,  it  is  assumed  that  all  the  activity  is  lofted.  For  buried  or  elevated 
bursts,  the  lofted  activity  is  decreased  by  the  following  fractions  derived  from 
dust  lofting  estimates  (ref.  4). 

F}  = 2.4  x 1 0~6 (l  1 2 . 5 + ((-9.11  x 10-12  H2  - 9.6  x 10-6)  H2 

+ 0.755)  h)2(32.7  + ((1.78  x 10-10  H2  - 2.52  x 10"5)  H2  + 0.851 ) h)  H < 0 

(34) 

Fz  = 8.6  x 10"8  (180-H)2  (360  + H),  0 < H < 180  (35) 

where 

H = 3.28  (HOB)  W */3 

HOB  = height  of  burst  (meters) 

The  activity  time  dependence  is 

f(t)  = 2.76  x 10  f1-23  -2.39  x 104  f1-45  (36) 

where  f(t)  is  normalized  to  one  at  3600  seconds. 

The  particle  size  distribution  and  the  activity  per  particle  of  size 
r microns  are  assumed  power  law  with  indices  n and  m respectively.  These 
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functions  are  combined  and  normalized  to  unity  between  a minimum  scale,  r • , 

and  a maximum  scale,  r „ 

max. 

H(r,  - <*♦  „ * »r^/  (Cx^  - (37) 

Nominal  values  for  m and  n are  2.7  and  -4  respectively.  Nominal  values  for 

r and  r . are  1000  and  0.5  microns  respectively.  The  particle  distribution 
max  min 

is  represented  by  a range  of  particles  of  decreasing  size.  The  activity  per 
micron  is  approximately 


A(r,t)  = 1500  f(t)  F M(r)  (roentgens  per  hour  per  KT  per  square  mile) 

(38) 

The  successive  particle  sizes  representing  the  particle  distribution  are 
chosen  starting  at  the  largest  size.  The  sizes  are  chosen  to  maintain  adequate 
resolution  over  the  particle  size  distribution  and  the  downwind  landing  positions 
of  particles  of  successively  smaller  size.  Let  us  define  a function,  T(r), 
which  is  approximately  proportioned  to  the  fall  time  of  a particle  of  size,  r, 
from  ten  kilometers. 


T(r) 


-1.9084 

= 6479.0  r 

-1  .3008 

= 383.1  r 


-0.6973 


= 6.1527  r 


, r < 105 
105  < r < 940 
, 940  < r 


(39) 


For  particle  indices,  j,  where 


r-i 


and 


nrj)  * 


Kp>  • ^.j)) 


NP  - j 


> B r. 

J 


(40a) 


r-1  [nr)]  = 


(40b) 


r^p  represents  the  largest  size  group,  and  NP  is  the  number  of  groups,  r. 


is  chosen  proportional  to  rj_-j 


rj  = rj-l/B 


(41) 
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Let  the  last  index  where  equation  (40)  is  true  be  m.  All  smaller  particle  group 
sizes  are  chosen  by 


rj 


T-i 


T(rm  ) * ( H 


(rM  - TQ-„)) 

NP  - m 


(42) 


The  parameters,  r-j , r^p,  NP,  and  8 are  chosen  to  maintain  the  specified  resolu- 
tion and  include  the  range  of  particle  sizes  important  in  a given  problem. 


The  particles  representing  the  particle  distribution  are  transported 
according  to  an  assumed  upward  fallout  cloud  velocity,  a calculated  fall-rate, 
and  a maximum  achievable  altitude  in  the  random  wind  field. 


Initially,  each  particle  is  assumed  to  rise  according  to  the  following 
equation. 

H(t)  = atb,  TST  < t < TSR  (43) 


where 

H(t)  = particle  altitude 

b = ALOG  (H./RST)/AL0G  (TSR/TST) 
a = (RST  + H0B)/TSTb 
RST  = initial  altitude 
TST  = initial  time 
TSR  = cloud  stabilization  time 
H.j  = maximum  achievable  altitude 

The  rise  velocity  of  a particle  is 

Vr(t)  = ab  tb"\  TST  < t < TSR 
= 0,  TSR  < t 


(44) 


The  particle  fall  rate  is  calculated  as  in  DELFIC5  with  an  extension  of  the  model 
for  large  particles.  Let 
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Co  = 

r (P  + 110.4)/(1 .458  x 10"6  pJ- 

5) 

ci  = 

c2  r cQ 

C3  = 

C0  C1  p 

(45) 

C4  = 

LOG10  (c3)  - 20.773 

P = 

air  density  (kg/m3) 

C2  = 

3.4  x 10-14 

The 

constant  c2  assumes  a 

fallout  density  of  2.6  g/cm3. 

The  fall 

rates  are 

Vf 

= (l  + 0.233/(rp))Cl  ( 

f 

41666.7  + c,  ('-2.3363  x 102 

k . 3 ' 

+ C3  (2.0154  - 

6.9105 

X 10'3  C3)))  * C3  < 140 

(46a) 

Vf 

= (l  + 0.233/(rp))  (50657.0  cj  c3(o.0011235(cu2  - 443 

.581), 

140  < c 

< 4.5 

x 107 

(46b) 

Vf 

= 0.294  (r/p)^ 

, c3 

> 4.5  x 10 

(46c) 

For  Vf  less  than  Vr,  the  particle  rises  with  the  rise  velocity.  As  the  particle 

rises,  V,  increases  and  V„  decreases.  After  V,  exceeds  V , the  particle  falls 
t r t r 

with  velocity  until  it  reaches  the  ground.  During  the  rise  and  fall  of  the 
particle,  it  is  transported  horizontally  according  to  the  random  winds  at  the 
particles  altitude  and  horizontal  position. 

The  maximum  achievable  altitudes,  , are  chosen  to  resolve  the  nuclear  cloud 
at  stabilization  altitude.  H.  varies  from  the  cloud  bottom,  ACB,  to  the  cloud 
top,  ACT.  The  number  of  intermediate  values  of  H..,  NLVL,  chosen  to  guarantee 
reasonable  downwind  resolution  of  the  fallout  pattern.  A particle  of  size  r is 
transported  for  each  value  of  H^.  To  account  for  the  activity  for  a given 
particle  size,  A(r,t)  of  equation  38  must  be  divided  according  to  H.  to  reflect 
the  division  of  the  fallout  cloud  into  layers.  In  addition,  it  is  necessary  to 
account  for  the  radius  of  the  cloud  which  is  needed  to  calculate  the  activity 
due  to  a given  particle  size  from  a given  layer  of  the  fallout  cloud.  The  cloud 
radius  at  time  t, , is 
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■.  ■[(«  4 - *tf 

where 

t = particle  land  tine 
t$  = time  when  V^  = Vr 
d = ALOG  (RMAX/RMIN)/ALOG(TSR/TST) 
c = RMAX/TSRd 

RMAX  = stabilization  radius 
RMIN  = radius  at  initial  time 
D = diffusion  constant 

The  activity  per  micron  given  for  t greater  than  t^  is 

A(r,t,Hi)  = 1500  f(t)  Fe  M(r)  G1/(nR|)  (48) 

where 

Gj  = (H2  - ACB)/(2(ACT-ACB)) 
gnlvl  = ' "mlwl-i )/(2(ACT-ACB)) 

Gi  = (Hi  + 1 -Hi  _ ^ ) / ( 2 ( ACT-ACB )) 

The  transport  now  can  be  summarized  by  the  following  steps.  A set  of 
particles  representing  the  particle  size  distribution  is  transported  with 
equal  to  ACB.  The  landing  points  define  a line  called  a trace.  Next,  another 
set  of  particles  are  transported  with  H.  equal  to  ACB  + DLVL  defining  a second 
trace.  DLVL  is  the  nominal  distance  between  successive  H^. 

DLVL  = (ACT  - ACB)/ (NLVL-1 ) (49) 

Next  for  each  landing  point  of  the  second  trace,  the  minimum  distance  to  the 
first  trace  divided  by  the  cloud  radius  is  calculated.  The  maximum  of  these 
values,  DMAX,  is  stored  in  an  array,  FLVL,  if  larger  than  one,  and  the  landing 
points  of  the  second  trace  are  permanently  retained.  If  DMAX  is  less  than  one 
the  last  trace  is  ignored.  In  either  case,  is  increased  by  DLVL  and  the 
next  trace  is  calculated.  Traces  are  calculated  for  the  selected  Hi  with  Hi 
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equal  to  ACT  being  the  last.  By  picking  reasonable  NLVL,  hence  DLVL,  this  pro- 
cedure results  in  a minimum  number  of  traces  representing  the  downwind  fallout 
pattern.  The  landing  points  defining  the  traces  are  stored  in  a transport  data 
base  for  each  burst  to  be  assessed  later  for  analysis. 

The  above  transport  procedure  can  become  expensive  if  there  are  a large 
number  of  bursts.  There  are  two  model  simpl ications  that  can  greatly  expedite 
the  calculations.  If  the  scenario  has  groups  of  similar  bursts  at  approximately 
the  same  time  within  a limited  area,  then  these  groups  can  be  represented  by  a 
single  group  cloud  whose  radius  is  equal  to  the  radius  of  the  area  of  the  burst 
group  plus  a single  burst  radius.  The  activity  of  the  new  cloud  is  increased 
according  to  the  number  of  bursts  in  the  group.  The  second  simpl ication  is  to 
use  an  abbreviated  transport  algorithm.  Three  traces  are  calculated  with  Hi 
equal  to  ACB,  (ACB  + ACT)/2,  and  ACT.  The  middle  traces  defines  the  landing 
positions  and  times  for  the  final  single  trace.  The  cloud  radius  for  each 
landing  point  on  the  middle  trace  is  the  average  of  the  distances  from  the 
first  and  third  trace  to  the  landing  point  plus  the  original  cloud  radius.  All 
of  the  activity  is  assigned  to  the  single  trace. 

Both  of  these  models  cost  some  accuracy  in  the  final  results.  The  group 
cloud  smears  out  the  fallout  in  the  local  area  of  the  bursts  and  the  loss  of 
resolution  is  determined  by  the  original  burst  locations.  Downwind,  the  group 
cloud  smearing  is  not  significant  compared  to  the  uncertainties  in  multiple 
cloud  dynamics.  The  simplified  transport  model  loses  resolution  within  the 
downwind  fallout  pattern.  This  model  does  deposit  approximately  the  correct 
amount  of  fallout  within  the  correct  downwind  area.  Thus,  the  simplified  model 
is  adequate  for  statistical  studies  and  limited  transport  times.  It  must  be 
abandoned  when  the  Increased  internal  resolution  of  specific  patterns  are  re- 
quired  or  the  transport  times  get  too  long.  The  resolution  available  with  the 

i 

multiple  level  transport  is  about  a local  cloud  radius.  This  is  the  best  avail- 
able with  any  disk  transport  model. 

In  the  preceding  discussion,  values  of  the  parameters  dealing  with  the  fall- 
out cloud  were  enumerated  but  not  numerically  defined.  The  fits  used  in  this 
model  are  as  follows: 
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TST  = 0.15  W0-333  (50a) 

RST  = 150.0  W0-333  + Maximum  (0.0,  HOB)  (50b) 

TSR  = 418.0  - 120.0  W0-5,  U < 1 

= 20.0(16.0  + 6.0  log  (W)  - W° • 2 7 ) , W > 1 (50c) 

10 

ACB  = 1800.0,  W < 1.5 

ACB  = 1447.0  wo.5385,  1.5  < W < 10.0  (50d) 

ACB  = 3263.0  W°-1854,  10.0  < W 

ACT  = 4674.0  W°-3080,  U < 10.0  (50e) 

ACT  « 6210.0  W°-1946,  W > 10.0 

NLVL  = (ACT-ACB)/1 000.0  (50f) 

RMAX  = 2553.0  W0-2310,  W < 30.0 

(50g) 

= 1884.0  W°-3203  W > 30.0 

RMIN  = 150.0  W0-333  (50h) 

D = 100.0  (50i) 
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SECTION  V 

FALLOUT  DATA  PROCESSING 

Fallout  data  is  usually  desired  in  one  of  several  formats.  Three  formats 
are  available  here.  The  first  format  involves  contours  of  the  mean  and  standard 
deviation  of  the  integrated  dose  at  specified  times  for  the  scenario  of  interest 
Contours  are  also  available  for  each  different  wind  pattern.  The  second  format 
supplies  the  man  and  standard  deviation  of  integrated  dose  for  selected  times 
for  an  inputted  list  of  targets.  The  last  format  supplies  the  statistics  of  the 
number  and  percentage  of  United  States  population  of  integrated  dose  at  geo- 
graphic points  supplied  by  a population  model  supplied  by  the  RAND  Corporation 
(ref.  6). 

All  three  formats  require  the  calculation  of  the  integrated  dose  at  a point 
from  which  trace  data  in  the  transport  data  base.  Let  us  consider  the  landing 
points  of  particles  of  sizes,  r^  and  r^,  at  times,  ti+^  and  ti , with  cloud 
radii,  Ri+1  and  Ri , with  activities,  A(ri+1,  ti+1 , H)  and  A (r.,  t^ , H),  on  a 
single  trace.  These  quantities  are  assumed  to  vary  linearly  between  the  end 
points  with  respect  to  time. 

Let  us  define  the  linear  dependence  as  follows: 


x 

y 

R 

A 


3 at  + b 
x x 

3 a t + b 

y y 

■ V * *>* 

■ <V  ♦ bA> 


f(t) 


where 


" (xi  ♦ i - xi>/(ti  + i - V 


bx  = Xi  - Vi 


(51) 
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ay  = (y,  * , - y,)/(t,  * , - t,l 


b = y . - a t . 
y y i 


aR  = (Ri  + 1 - "i^i  + 1 - V 
bR  = Ri  - aRt. 

aA  = A^ri  + 1»  ti  + 1*  + ])  - A(ri»  V H)/f(ti)/(ti  + 1 - tj) 

bA  = A(ri,  tr  H)/f(ti)  - aAti 

In  order  for  a point  (xQ,  yQ)  to  be  affected  by  fallout,  the  following 
condition  must  be  true  for  some  t between  t^  and  t^  + 

(x  - x0)2  + (y  - yQ)2  < R2 

*2  [ai  * a)  - aft]  * l[2ax  (Bx  - xo>  * 2ay  (by  x V - 2aRbA] 


[' 


bx  ' V2  * (by  * y 


o>2  - Hi]  ■ 


(52a) 


(52b) 


In  general,  equation  (52b)  has  zero,  one,  or  two  roots  for  t.  If  there  is 
one  or  less  roots,  then  no  fallout  has  landed  at  (xQ,  yQ).  If  there  are  two 
roots,  then  the  fallout  depends  on  the  values  of  the  roots.  Let  the  large  and 
small  roots  be  t£  and  t$  respectively.  Then  redefine  these  quantities. 


t^  = minimum  (t^,  ti  + 1 , T) 
ts  = maximum  (t$,  t- ) 


(53) 


where  T is  the  sample  time  for  the  integrated  dose.  If  t$  is  greater  or  equal 
to  t?,  then  no  fallout  occurs.  If  not,  then  there  is  some  contribution  to  the 
fallout.  Assuming  that  the  fallout  is  uniformly  distributed  in  the  horizontal 
direction,  the  contribution  to  the  integrated  dose  at  xQ,  yQ  is 


• (m)  /-  (■  / 


/ du  u f (u)  + bA  J du  f (u)  1 (54) 
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The  first  factor  on  the  right  side  is  used  to  maintain  proper  normalization 
because  the  outer  integral  is  better  thought  of  as  an  integral  over  particle 
size.  Ds  is  calculated  for  each  successive  pair  of  locations  along  each  trace 
to  get  the  total  integrated  dose. 

To  maintain  proper  resolution  transverse  to  the  general  direction  of  the 
trace,  the  traces  must  not  be  separated  by  more  than  the  local  cloud  radius. 

The  traces  in  the  transport  data  base  are  not  that  close.  To  maintain  adequate 
resolution,  pseudo  traces  are  linearly  interpolated  between  the  calculated 
traces.  This  procedure  maintains  reasonable  resolution  without  having  to 
transport  an  excessive  number  of  particles.  Where  the  simplified  transport 
is  used,  no  interpolation  is  required  as  there  is  only  one  trace. 
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SECTION  VI 
RESULTS 

Two  scenarios  were  run  to  demonstrate  the  effects  of  wind  variability.  The 
first  scenario  involved  7000  one-megaton  warheads  targeted  on  a proposed  advanced 
missile  (MX)  complex  near  Yuma,  Arizona.  The  second  scenario  targeted  the 
Minuteman  Missile  (MM)  complexes  on  the  great  p 1 <3 i n s with  approximately  two 
one-megaton  weapons  per  missile  silo.  In  both  cases  group  clouds  were  used  for 
each  local  complex.  Thus  one  cloud  modeled  the  MX  problem  and  six  clouds,  one 
for  each  wing,  were  used  for  MM.  The  cloud  radii  were  all  set  at  sixty  miles. 

The  remaining  cloud  dynamics  were  identical  to  a single  one-megaton  burst  at 
ground  level.  The  fission  fraction  of  all  bursts  was  0.5.  The  particle  size 
distribution  and  other  necessary  details  were  discussed  earlier  in  the  trans- 
port or  wind  sections.  The  wind  statistics  used  were  for  January.  The  correla- 
tion lengths  and  other  statistical  properties  of  the  wind  field  were  previously 
discussed.  Figure  1 shows  the  mean  winds  for  January  at  two  kilometers  altitude. 
The  winds  trend  toward  the  east  as  expected.  Figure  2 shows  a sample  wind 
pattern  at  the  same  altitude  where  the  variability  effects  have  been  added. 

This  pattern  shows  typical  circulatory  wind  patterns  characteristic  of  real 
weather  patterns.  It  is  the  effects  of  these  circulatory  features  that  are 
under  study. 

Figures  3-7  show  sample  dose  patterns  for  five  different  wind  patterns  for 
the  MX  scenario.  Figure  8 shows  the  fallout  pattern  using  the  mean  winds. 
Comparison  shows  that  real  winds  can  carry  fallout  to  regions  that  would  not  be 
covered  by  mean  winds.  Figure  9 shows  the  average  dose  contours  for  the  two 
hundred  cases  run  for  the  MX  scenario.  The  roughness  of  the  contours  indicate 
that  even  the  large  number  of  cases  used  give  only  moderate  averages  at  the 
low  dose  levels  indicating  the  residual  effects  of  high  variability.  A more 
useful  comparison  between  mean  and  variable  winds  is  in  table  1.  The  table 
shows  the  number  and  percent  of  the  United  States  population  exposed  to  the 
listed  dose  for  the  mean  and  variable  winds  with  the  standard  deviation  for  the 
variable  winds.  The  most  obvious  significantly  underestimated  by  the  mean  wind 
calculation.  This  is  caused  by  the  winds  taking  the  patterns  occasionally  over 
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heavily  populated  areas  raising  the  mean  expected  doses.  This  effect  is  clearly 
accentuated  by  the  high  strength  fallout  source  in  a localized  area  which  leads 
to  very  long  intense  fallout  patterns.  Figure  10  demonstrates  the  most  extreme 
case  out  of  the  two  hundred  runs.  This  case  exposed  over  fifty  million  people 
to  more  than  1300  R.  This  very  high  casualty  count  was  caused  by  three  factors. 

First,  the  fallout  cloud  did  not  appreciably  disperse  laterally.  Second,  the 
cloud  rapidly  moved  to  the  east.  Lastly,  the  pattern  went  over  the  heavily 
populated  eastern  corridor.  This  is  a very  extreme  case  as  shown  by  fig.  11 
which  shows  the  cumulative  casualty  distribution  for  the  MX  scenario.  There  is 
only  about  a three  percent  chance  that  the  population  exposed  to  over  1300  R 
would  exceed  ten  million.  The  skewness  of  the  cumulative  distribution  results 
in  the  mean  casualties  at  a given  level  as  shown  in  table  1 differing  from  the 
50  percent  probability  casualty  count  in  figure  11.  Also  shown  is  a comparison 
with  a RAND  calculation  (ref.  6)  at  the  1300  R level.  The  agreement  is  good 
considering  the  modeling  uncertainties. 

Figures  12-16  shows  sample  dose  patterns  for  the  MM  scenario.  The  patterns 
from  each  field  are  now  much  shorter  because  the  sources  are  now  spread  out  with 
less  total  yield.  Figure  17  shows  t.  mean  wind  calculation.  Again  the  varia- 
bility causes  regions  to  be  exposed  that  otherwise  would  not  be.  Figure  18 
shows  the  mean  dose  averaged  over  fifty  different  winds.  The  smooth  contours 
with  less  cases  clearly  demonstrates  the  decreased  impact  of  the  wind  variability 
with  dispersed  sources.  Table  2 repeats  table  1 for  the  MM  calculations.  In 
this  case  comparison  between  the  mean  wind  dose  and  the  mean  dose  over  winds  show 
little  difference.  Figure  19  shows  the  cumulative  casualty  distribution  for  MM. 

The  curves  are  nearly  linear  over  most  of  the  range  which  is  consistent  with  the 
expected  casualties  at  a given  level  on  table  2,  being  approximately  equal  to 
the  50  percent  probability  casualty  count. 

The  two  scenarios  indicate  that  the  impact  of  variability  decreases  as  the 
fallout  sources  are  dispersed.  For  the  MX  scenario,  the  intense  concentrated 
source  resulted  in  long  patterns  which  lead  to  high  variability  in  the  casualties 
as  the  patterns  wandered  over  large  areas  of  the  country.  For  the  more  dispersed 
MM  scenario,  the  individual  patterns  were  smaller,  decreasing  the  variability. 

The  implication  is  clear  that  for  many  widely  dispersed  sources  such  as  massive 
laydowns,  the  wind  variability  effects  can  be  ignored  and  mean  winds  used  for 
casualty  assessments.  Only  in  localized  or  very  intense  scenarios  with  long 
patterns  does  wind  variability  appear  to  be  significant. 
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Figure  1.  Mean  January  ,;'nds 
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Figure  7.  Sample  MX  Fallout  Dose  Contours  at  T = 108 


29 


AFWL-TR-78-64 


CONTOUR  SCALE 
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Figure  8.  MX  Mean  Wind  Dose  at  T = 108 
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CONTOUR  SCALE 
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Figure  M.  Sample  MM  Fallout  Dose  Contours  at  T = 108 
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CONTOUR  SCALE 
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Figure  16.  Sample  MM  Fallout  Dose  Contours  at  T = 108 


AFWL-TR-78-64 


1 1 . 50E*02 

2 4.50E+02 

3 1.35E+03 


NUMBER/PERCENT  OF  POPULATION  EXPOSED  AT  TIME  = 108 


AFWL-TR-78 


1000000 


CASUALTIES  (MILLIONS) 


AFWL-TR- 


20 


10 


0.0  0.2  0.4  0.6  0.8  1.0 
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Figure  19.  MM  Cumulative  Casualtv  distribution 
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